# Copyright 2013 DEVSIM LLC
#
# SPDX-License-Identifier: Apache-2.0

from devsim import (
    add_2d_contact,
    add_2d_mesh_line,
    add_2d_region,
    contact_equation,
    contact_node_model,
    create_2d_mesh,
    create_device,
    edge_from_node_model,
    edge_model,
    element_from_edge_model,
    equation,
    finalize_mesh,
    get_contact_charge,
    node_model,
    node_solution,
    set_parameter,
    solve,
    write_devices,
)

device = "MyDevice"
region = "MyRegion"

xmin = -25
x1 = -24.975
x2 = -2
x3 = 2
x4 = 24.975
xmax = 25.0

ymin = 0.0
y1 = 0.1
y2 = 0.2
y3 = 0.8
y4 = 0.9
ymax = 50.0

create_2d_mesh(mesh=device)
add_2d_mesh_line(mesh=device, dir="y", pos=ymin, ps=0.1)
add_2d_mesh_line(mesh=device, dir="y", pos=y1, ps=0.1)
add_2d_mesh_line(mesh=device, dir="y", pos=y2, ps=0.1)
add_2d_mesh_line(mesh=device, dir="y", pos=y3, ps=0.1)
add_2d_mesh_line(mesh=device, dir="y", pos=y4, ps=0.1)
add_2d_mesh_line(mesh=device, dir="y", pos=ymax, ps=5.0)

device = device
region = "air"

add_2d_mesh_line(mesh=device, dir="x", pos=xmin, ps=5)
add_2d_mesh_line(mesh=device, dir="x", pos=x1, ps=2)
add_2d_mesh_line(mesh=device, dir="x", pos=x2, ps=0.05)
add_2d_mesh_line(mesh=device, dir="x", pos=x3, ps=0.05)
add_2d_mesh_line(mesh=device, dir="x", pos=x4, ps=2)
add_2d_mesh_line(mesh=device, dir="x", pos=xmax, ps=5)

add_2d_region(
    mesh=device, material="gas", region="air", yl=ymin, yh=ymax, xl=xmin, xh=xmax
)
add_2d_region(mesh=device, material="metal", region="m1", yl=y1, yh=y2, xl=x1, xh=x4)
add_2d_region(mesh=device, material="metal", region="m2", yl=y3, yh=y4, xl=x2, xh=x3)

# must be air since contacts don't have any equations
add_2d_contact(
    mesh=device, name="bot", region="air", yl=y1, yh=y2, xl=x1, xh=x4, material="metal"
)
add_2d_contact(
    mesh=device, name="top", region="air", yl=y3, yh=y4, xl=x2, xh=x3, material="metal"
)
finalize_mesh(mesh=device)
create_device(mesh=device, device=device)


###
### Set parameters on the region
###
set_parameter(device=device, region=region, name="Permittivity", value=3.9 * 8.85e-14)

###
### Create the Potential solution variable
###
node_solution(device=device, region=region, name="Potential")

###
### Creates the Potential@n0 and Potential@n1 edge model
###
edge_from_node_model(device=device, region=region, node_model="Potential")

###
### Electric field on each edge, as well as its derivatives with respect to
### the potential at each node
###
edge_model(
    device=device,
    region=region,
    name="ElectricField",
    equation="(Potential@n0 - Potential@n1)*EdgeInverseLength",
)

edge_model(
    device=device,
    region=region,
    name="ElectricField:Potential@n0",
    equation="EdgeInverseLength",
)

edge_model(
    device=device,
    region=region,
    name="ElectricField:Potential@n1",
    equation="-EdgeInverseLength",
)

###
### Model the D Field
###
edge_model(
    device=device, region=region, name="DField", equation="Permittivity*ElectricField"
)

edge_model(
    device=device,
    region=region,
    name="DField:Potential@n0",
    equation="diff(Permittivity*ElectricField, Potential@n0)",
)

edge_model(
    device=device,
    region=region,
    name="DField:Potential@n1",
    equation="-DField:Potential@n0",
)

###
### Create the bulk equation
###
equation(
    device=device,
    region=region,
    name="PotentialEquation",
    variable_name="Potential",
    edge_model="DField",
    variable_update="default",
)


###
### Contact models and equations
###
for c in ("top", "bot"):
    contact_node_model(
        device=device, contact=c, name="%s_bc" % c, equation="Potential - %s_bias" % c
    )

    contact_node_model(
        device=device, contact=c, name="%s_bc:Potential" % c, equation="1"
    )

    contact_equation(
        device=device,
        contact=c,
        name="PotentialEquation",
        node_model="%s_bc" % c,
        edge_charge_model="DField",
    )

###
### Set the contact
###
set_parameter(device=device, name="top_bias", value=1.0e-0)
set_parameter(device=device, name="bot_bias", value=0.0)


edge_model(device=device, region="m1", name="ElectricField", equation="0")
edge_model(device=device, region="m2", name="ElectricField", equation="0")
node_model(device=device, region="m1", name="Potential", equation="bot_bias;")
node_model(device=device, region="m2", name="Potential", equation="top_bias;")


# solve -type dc -absolute_error 1.0 -relative_error 1e-10 -maximum_iterations 100 -solver_type iterative
solve(
    type="dc",
    absolute_error=1.0,
    relative_error=1e-10,
    maximum_iterations=30,
    solver_type="direct",
)

element_from_edge_model(edge_model="ElectricField", device=device, region=region)
print(get_contact_charge(device=device, contact="top", equation="PotentialEquation"))
print(get_contact_charge(device=device, contact="bot", equation="PotentialEquation"))

write_devices(file="cap2d.msh", type="devsim")
write_devices(file="cap2d.dat", type="tecplot")
write_devices(file="cap2d", type="vtk")
